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Analysis of discretization errors in LES 

By Sandip Ghosal 1 

1. Motivation and objectives 

All numerical simulations of turbulence (DNS or LES) involve some discretization 
errors. The integrity of such simulations therefore depend on our ability to quantify 
and control such errors. In the classical literature (see e.g. Chu 1978) on analysis of 
errors in partial differential equations, one typically studies simple linear equations 
(such as the wave equation or Laplace’s equation). The qualitative insight gained 
from studying such simple situations is then used to design numerical methods 
for more complex problems such as the Navier-Stokes equations. Though such an 
approach may seem reasonable as a first approximation, it should be recognized 
that strongly nonlinear problems, such as turbulence, have a feature that is absent 
in linear problems. This feature is the simultaneous presence of a continuum of 
space and time scales. Thus, in an analysis of errors in the one dimensional wave 
equation, one may, without loss of generality, rescale the equations so that the 
dependent variable is always of order unity. This is not possible in the turbulence 
problem since the amplitudes of the Fourier modes of the velocity field have a 
continuous distribution. The objective of the present research is to provide some 
quantitative measures of numerical errors in such situations. Though the focus of 
this work is LES, the methods introduced here can be just as easily applied to DNS. 
Errors due to discretization of the time- variable are neglected for the purpose of 
this analysis. 

2. Accomplishments 

In this report, analytical expressions for the power spectra of errors due to the 
spatial discretization of the Navier-Stokes equations are derived. In § 2.1, an ex- 
pression for the numerical error is presented as the sum of “finite-differencing”, 
“aliasing”, and “modeling” errors that have different origins. In § 2.2, expressions 
for the power spectra of the first two kinds of errors are derived as well as the corre- 
sponding expressions for the subgrid and total nonlinear terms. The essential tool 
that makes the derivation of such an analytical expression possible is the “joint- 
normal hypothesis” for turbulent velocities. The essential technique is identical to 
that used by Batchelor in his derivation of the pressure spectrum of turbulence 
from the energy spectrum (Batchelor 1951, 1953). These results are applied to the 
LES of turbulence in § 2.3 to obtain some measure of numerical errors in finite- 
difference schemes, which are increasingly being used in turbulence computations 
on flows with complex boundaries. This report summarizes the essential results, 
the details of the mathematical development will be presented elsewhere (Ghosal 
1995 - henceforth referred to as “paper 1”). 

1 Present address: CNLS (MS-B258), LANL, Los Alamos, NM 87545 
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2.1 Calculation of discretization errors 

Any representation of the true velocity field in a turbulent flow on a finite grid 
is necessarily approximate. One must be careful to distinguish between errors due 
to the finiteness of the representation and the “discretization error” of a numeri- 
cal scheme. In a numerical simulation, the velocity field at any time-step can be 
regarded as an element of a vector space with a finite number of dimensions (N) 
where N is the number of variables retained in the computation. This is an ap- 
proximate representation in a subspace of the larger vector space that contains the 
true solution. The best possible approximation to the true solution in the subspace 
is the projection onto the subspace (in fact that is the definition of a projection 
operator — see e.g. Helmberg 1969). The “ideal” or “best approximation” to the 
Navier-Stokes operator in the finite subspace is that operator that ensures that the 
numerical solution remains “locked” to the projection of the true solution at all 
times as both vectors move around in their respective vector spaces. It may be 
shown (see paper 1) that this condition is achieved by spectral methods (or prop- 
erly dealiased pseudo- spectral methods) in the absence of subgrid modeling errors. 
By “discretization error”, E, of a numerical method we mean the deviation of the 
right hand side evaluated with the method from what would have been obtained 
if the right-hand side of the full Navier-Stokes equation were projected into the 
computational subspace. Thus, for a spectral method used in conjunction with an 
‘exact’ subgrid model, E = 0. 

In order to evaluate the formal expression for the error E, one needs to intro- 
duce a basis. The most advantageous choice is the 3D Fourier-basis since in Fourier 
space differentiations reduce to multiplication by wavevectors and numerical differ- 
entiation reduces to multiplication by modified wave vectors (see e.g. Vichnevetsky 
1982). We now restrict our attention to flows in a periodic cubical box. Further, 
while considering finite-difference methods, the grid will be assumed uniform in ev- 
ery direction. Let Ei( k) denote the components of E in the Fourier-basis with i = 1 
to 3 corresponding to the x, y, and z directions respectively. Then £i(k) can be 
written as 

£,(k) = £< FD) ( k) + £< alias) (k) + E^ modei \k). (1) 

The first term arises because of the inability of the finite-differencing operator, 
6/8x *, to accurately compute the gradient of short- wavelength waves. We call this 
the “finite-differencing error.” It vanishes for a spectral method that can differenti- 
ate waves of all wavelengths exactly. The second term arises due to the method of 
computation of the nonlinear term by talcing products in physical space on a dis- 
crete lattice. This is called the “aliasing error” and is well known in the literature 
on pseudo-spectral methods (Canuto et al. 1988, Rogallo 1981). The last term is 
the difference between the true subgrid force and that computed using a subgrid 
model. We call this the “modeling error” . 

In the following analysis, the magnitude of the error E will be characterized by 
statistical properties such as its power spectral density. Such statistical measures 
can be precisely defined only in the limit where the wavevector can assume a con- 
tinuum rather than a discrete set of values. In physical space this implies that we 
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are considering the grid size A and some characteristic scale of turbulence A fixed 
and taking the limit as the size of the box L — > oo. In actual simulations, of course, 
the box size L is finite. However, L is taken much larger that A or A so that smooth 
power spectra can be defined and computed statistical quantities are not changed 
when the box size is increased further. This ensures that the computed quantities 
are indistinguishable from the ideal limit, L — ► oo. For the purpose of theoretical 
analysis it is advantageous to take the limit L —► oo first rather than at the end of 
the computation. Thus, in the Fourier-basis, the exact solution will be characterized 
by a continuum of wave vectors k <E R 3 and the numerical solution will be character- 
ized by the subset k G B where B = [-Jfc™ ax , Jfc™ ax ] x [-fc™ ax , Jb™ ax ] x [-jfe™ 1 ax , k?* x \. 
We will assume for simplicity that the grid length A is the same in all three di- 
rections so that A:™ ax = A:™ ax = fc™ ax = k m = 7r/A. Further, we will consider 
the LES “filter-width” and the grid length to be identical. This condition will be 
relaxed in § 2.3.3. In the limit of infinite box size, the discrete Fourier transform 
and its inverse take the form (a factor of L z / 87r 3 is ‘absorbed’ in the definition of 
the transform) 

<K k ) = E ex P(-* k • x ) ^( x ) = J b <*M k ) exp(i' k • x) (2) 

where the summation is over all lattice points over the infinite cubic lattice of 
spacing A and the integration over wave space ranges over all vectors k € B. The 
following useful identity is readily derived by taking the limit of infinite box size: 

E ex p< iK • x ) = E - a ) ( 3 ) 

X a€A 

where is the Dirac delta function, A is the set of wavevectors of the form 
(2p& m , 2gfc m , 2rk m ) where p, q and r are integers (positive, negative or zero), and K 
is any vector (not necessarily restricted to B). [This relation is familiar in solid state 
physics (see e.g. chapter 1, pg. 12 of Jones & March 1973) where the set A goes by 
the name “Reciprocal Lattice”.] When the lattice spacing A — ► 0, the summation 
over lattice points in (2) becomes an integral over space and the usual continuous 
Fourier-transform is recovered. In this limit, the right hand side of relation (3) be- 
comes simply £(K) and (3) reduces to the familiar expansion of the delta-function 
in terms of exponentials. 

Let us first consider the effect of projecting the exact right-hand side of the Navier- 
Stokes equation onto the Fourier-basis with wavevector k. The ‘zth component’ is 
given by 

-*P,mn( k ) [ jT J B dk'dk"6(k' + k " - k )« m ( k >n( k ") + W k ) - ^ 2 «,( k ) 

where Pi mn has its usual meaning (see e.g. Lesieur 1987) and r mn is the exact 
subgrid stress in Fourier-space. The Einstein summation convention for tensor 
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indices is implied throughout this report except where otherwise noted. If the 
exact derivative operator d/dxk is replaced by the numerical differentiation 6/Sxk , 
multiplication by wavevectors k are replaced by multiplication by the corresponding 
modified wavevectors k. Thus, we obtain 


£, (FD) (k) =* [Pimni k) - Pimn{ k)] J f dk'dk W tf(k' + k" - k)u m (k')« B (k") 


+ i 


P imn (k)f" n (k) - P im „(k)f" w (k)] + p(k 2 - * 2 )u,(k), 


( 4 ) 


where T*f n ( k) is the “modified subgrid model” obtained by replacing all multipli- 
cation by wavevectors (if any) in the subgrid model f^ n (k) by the corresponding 
multiplication by modified wavevectors. 

To obtain the aliasing error, we consider the effect of evaluating the nonlinear 
term in physical space: 




u m (x)u n (x) + T £* n (k) 


On using the definition (2) of the discrete Fourier transform we have 
u m (x)u„(x) = ^2 u m(x)«n(x)exp(— ik • x). 

X 

When u m (x) and u„(x) in (5) are expanded in the Fourier-basis we get 
CftKW = £^Y,f B f B rfk ' dk "«m(k')ti n (k")exp[i(k' + k" - k) • x]. 


( 5 ) 


( 6 ) 


The summation over lattice points can be performed using (3), 

f dk'dk"u m (k> n (k>(k' + k" - k - a). (7) 

a€A Jb 


All the terms in the sum over a £ A with the exception of a = 0 are clearly “spurious 
contributions” and constitute the ‘aliasing error’. Thus, we have 

E^ 9 \k) = tP imn (k)T f f dk'dk"6(k' + k" — k — a)u m (k')u„(k ,/ ) + 6f„ n (k) 

»eAo jBjB 

( 8 ) 

where Ao consists of the vectors (2pfc m , 2qk m , 2rk m ) where p, q and r can indepen- 
dently take on the values 0 or ±1 but excluding the case p = q = r = 0. The reason 
integer values of p, q and r with modulus greater than 1 are not included in Ao is 
that the relation a = k' + k" — k cannot be satisfied for such values if k, k', k" £ B 
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and hence the delta function ensures that they do not contribute to the sum. The 
last term ST^ n ( k) is the contribution to the aliasing error from the subgrid model. 
Obviously it depends on the nature of the model. For a subgrid model that uses 
a constant eddy viscosity, this term is linear in the resolved fields, and hence there 
is no contribution to the aliasing error. For more complicated models such as the 
Smagorinsky model, it is difficult to evaluate this contribution analytically due to 
the complicated nature of the nonlinearity. 

The expressions (4) and (8) for the finite-differencing and aliasing errors involve 
the subgrid model t-j . Modeling errors associated with subgrid models are difficult 
to estimate, and, further, there is no obvious way to single out for this study any 
one among the wide variety of subgrid models in use. It is therefore advantageous to 
separate the issue of subgrid modeling from the issue of discretization errors which 
is the subject of this paper. In order to accomplish this, we introduce the concept 
of the “ideal subgrid model”: 

T ij=Tij(x,t) (9) 

where 7y/(x,t) is the exact subgrid stress. One might think of the “ideal subgrid 
model” (9) in the following way. Imagine a DNS with an infinitely greater resolution 
running concurrently with the given LES. At every time-step the exact subgrid stress 
is computed from the DNS fields and supplied to the LES simulation as a function 
of position. The rest of the analysis in this paper will be presented for such an 
idealized LES. Since tJj is already given as a function of position and time and 
involves no computation, it does not contribute to aliasing errors. Thus, for such 
an idealized LES, ST^ n = 0 in (8). The contribution from the subgrid model is not, 
however, zero for the finite-differencing errors even for the ideal model (9). This is 
because the model is (inaccurately) differentiated for computing the pressure and 
the subgrid force. The subgrid terms in (4) for the ideal model (9) are given by 

_ rpM _ - Thus 
7 mn — x mn ~ T mn • J-IlUb, 


E 


(FD) 


(k) =* [p, mn (k) - P im „(k)] J j dk'dk"<5(k' + k" - k)tt m (k')fi„(k") 

+ v(k 2 - 


( 10 ) 


The integration in (10) now ranges over the entire wave-space. Clearly, for this 
ideal model 

= *P,mn(k) [4">(k) - f mn (k)] = 0. (11) 

2.2 Power spectra 

In this section, analytical expressions for the power spectra of the finite-differencing 
error, aliasing error, subgrid and total nonlinear term are derived. 

2.2.1 Finite- differencing error 

The power spectrum of the finite- differencing error is defined by £( FD \k), where 


£< FD >(fe) 

47 rk 2 




( 12 ) 
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{ } Q denotes angular average in wave-number space over the surface of the sphere 
|k| = k and V is the volume of the physical box containing the fluid. 

From (4), we have, 


(£f D (k)£f D (k)*) = 

A, mn (k,k)A* p? (k,k) J dk'dk"(u m (k')u n (k - k>;(k>;(k _ k ")) 


+ 2 


*A* mn (k, k)(P - k 2 ) j dk'(^(k'K(k-k')Mk)) 


(13) 


+ I , 2 |P-fc 2 | 2 (ti i (k)u*(k)) 


where { ) denotes ensemble average, * denotes complex conjugate, 3 denotes the 

imaginary part, and A; mn (k,k) = Pi mn (k)-Pj mn (k). The following two properties 
of the A imn tensors follow immediately from the corresponding properties of Pimn ; 

A imm — 0? A tmn = A inm . 

In order to make further analytical work possible with (13), we now introduce 
the “Millionshchikov hypothesis” (see e.g. Monin and Yaglom 1979) that in fully 
developed turbulence, the joint probability density function of any set of velocity 
components at arbitrary space- time points can be assumed to be joint-normal. The 
joint-normal hypothesis was originally evoked in turbulence in an attempt to close 
the hierarchy of equations for moments (see e.g. Lesieur 1987). Though this did 
not succeed, the joint-normal hypothesis has been successfully used in other con- 
texts. Thus, Batchelor (Batchelor 1951) used it with success to predict the pressure 
spectrum of isotropic turbulence. The joint-normal hypothesis implies in particular 

(tZj(Xi)uj(x 2 )u*(x 3 )ti/(x 4 )) =(u,(Xi)Uj(x 2 ))(u*(X 3 )u,(x 4 )) 

+(u,(xi)u*(x 3 ))(ti>(x 2 )u/(x 4 )) (14) 

+(«.(x i)u/(x 4 ))(u>(x 2 )u*(x 3 )) 

and that all third order moments are zero. Here u(x,f) is the true velocity field 
defined at all space time points. On talking the (continuous) Fourier transform of 
(14) and assuming the turbulence to be homogeneous, we have, 

(u,(ki)u J (k 2 )ujt(k 3 )u/(k 4 )) =<5(ki + k 2 )£(k 3 + k 4 )$ij(k 2 )$fc/(k 4 ) 

+tf(k, + k 3 )£(k 2 + k 4 )<Mk 3 )<Mk 4 ) (15) 

+tf(k 4 +k 4 )*(k 2 +k 3 )S i i(k 4 )* i *(k 3 ), 

where^ij is the Fourier transform of the correlation tensor i? tJ (x 2 — Xi) = (w,(xi )u ; (x 2 )). 
On substituting (15) into (13) we get after some algebra (see paper 1) 

€ (FD) (k) = |8»lb 2 A imB (k,k)A jM , (k,k) J ^ p (k')$* n ,(k-k')d 3 k' 

+v 2 \k 2 -k 2 \ 2 * ii (k)} Q . 


(16) 
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Eq. (16) is the general result for homogeneous turbulence. If in addition, the 
turbulence is isotropic, simplifies (Batchelor 1953) to 

= <17) 

where E(k ) is the three dimensional energy spectrum and 6ij is the Kronecker-delta 
symbol. The integral in the first term of (16) may be written after substitution of 
(17) as 


•W k) sSrf 2 J *;,(k')*;,(k - k') 

e f E(P)E(Q) . 2 2 , 

— 2tt y P 4 Q 4 ^ ^ °mpOnq “ FmFpQ °nq ~ 

x6(P + Q-k) d 3 Pd 3 Q. 

This integral can be simplified (see paper 1). The result is 


QnQqP 2 fi mp + PmPpQnQ, 


(18) 


Jmpnq(k) — -^1 (^)^mp^ng H~ 7*2 (^)(^mn^pj ”h &pn&mq^ 


+ F 3 (k) 


kmkp - 


nq 


+ 


knkq c 
~fcT dm P 


+ F 4 (k) 


where 


F 1 (k) = ^ [7I 4 + 6/ 3 - 2I 2 + 5 /,] 

lo 

F 2 (k) = ^ [-3/4 + 2/3 + 2 / 2 - /,] 
Fs(k) = ^ [-15/4 - 6/3 + 2 / 2 + 3A] 
Fiik) = ^ [45/4 - 30/3 - 6/2 + 7/,] . 


kmkpk n kq ( 19 ) 

k 4 


( 20 ) 


The terms I m are defined as 


/m — k 



d V E(kOE(kr,)W m ((,ri) 


where the weights W m axe defined as follows: 


Wi(t,v) 

W 2 (Z,r,) 

w 3 u,v) 

W 4 (t, V ) 


1_ 

tv 

(1 -C 2 -*? 2 ) 2 

4£ 3 j/ 3 

a ±e_zvy 

^V 

[l - (£ 2 - >? 2 ) 2 ] 2 

16t 3 V 3 


( 21 ) 


( 22 ) 
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Therefore, after substituting (19) in (16) and using the properties A i mm = 0 and 
A* mn = A, nm , the following expression is obtained for the power spectrum of the 
finite- differencing error (no summation over repeated indices!): 


£< FD >(fc) 

[F,(*) + F 2 (*)J { E I Aj mn (k,k)|^ +2 F 3 (k){ 

Q 


i ,m,n 


^,2 ^*"»n(k, k)A ipn (k, k) 

i,m,n,p 


+ F 4 (k){ 


E 


— A <mn(k , k)A* pg (k, k) \ + 2u 2 E(k){\k 2 - k 2 ' 2 


A : 4 




n 


(23) 

In (23), the functions Fi(k), F 2 (k ), Fs(k) and F^k) are known once the energy 
spectrum is specified. They are not affected by the choice of numerical schemes. 
On the other hand, the coefficients of these functions in (23) depend only on the 
numerical method (through the dependence of k on k) and are quite independent 
of the physical spectrum. Thus, given a specific numerical scheme and energy 
spectrum, (23) can be used to compute the power spectrum of the finite-differencing 
error. This is done in §2.3 for various representative numerical schemes. 

2.2.2 Aliasing errors 

The power spectrum of the aliasing error is defined by 


£( alias >(jfc) 

47rA; 2 

From (8) one obtains 


.. 87 r 3 

= lim — 
v —+00 V 


{(£< alias) (k)£j^ ias) (k)*)} o . 


(24) 


(Ef Us (k)F^(k) m ) = P imn (k)P; p9 (k) Y, f f f f dk,dk 2 dk 3 dk 

a,a'GA 0 ^ 8 ^ B *' B B 


(25) 


x (u m (ki)u n (k2)up(k3)u!(k 4 ))<5(k + a-ki - k 2 )6(k + a' - k 3 -k 4 ). 


On applying the joint-normal hypothesis, (15), one gets after some algebra (see 
paper 1) 

£ (alias) (Jfc) = 

8nk 2 Y {Pimn(k)P; pq (k) f j dk'<flc"$^p(k')$* g (k")6(k + a — k' — k")l . 
a€A 0 1 JbJb > « 

(26) 

The integral in (26) is difficult to handle analytically because integration over the 
cubical region B destroys the spherical symmetry of the problem exploited in the 
computation of S FD (k) in the last section. In order to make analytical progress, the 
following approximation is introduced. The region which is a cube in k-space, 
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is replaced by the largest sphere contained in it. Clearly, this procedure can be 
implemented simply by removing the suffix ‘5’ from the integral signs in (26) and 
replacing the energy spectrum E(k ) by 

jfc) = | if ^ < k m fo7\ 


otherwise. 


The superscript ‘min’ indicates that this procedure underestimates the true aliasing 
error by failing to take account of the contribution of modes close to the eight corners 
of the cube. An alternative method that overestimates the error can be provided by 
replacing the cube by the smallest sphere that contains it. To obtain this estimate 
one needs to use in place of E miTl the following spectrum; 


__ f E(k) if k < y/3 k n 
l 0 otherwise. 


The true aliasing error is then expected to lie between these two bounds. With 
the approximation so described, and with the energy spectrum defined as in (27) 
or (28), the integral in (26) may be extended to the entire wave space. Thus, one 
obtains 

£^™\k) = £ {p, mn (k)P* ? (k) J mpni (k + a)}^ . (29) 

a6Ao 

Substitution of the expression for <7 mpn? gives (no summation over repeated in- 
dices!): 


£<**->(*)= £ \[F 1 (K) + F 2 (K)} J2 |P,mn(k )| 2 


K m K, 


+2F 3 (K) £ — ^j- E -Pi mn (k)P)p n( k) -1- F 4 (K) £ 


K m K p K n K, 


P imn (k)P*,(k) 


where K = k + a. Note that in this case the Fi(K ) does depend on the direction 
of k so that the Fi(K) cannot be extracted from the { }q operation. Though the 
summation over the set Ao consists of 3 3 — 1 = 26 terms, for a cubical box one only 
needs to evaluate 3 terms due to symmetry. Indeed, the full set of “aliasing modes”, 
a € Ao, fall into three classes (Rogallo 1981): 

'(±2fc m ,±2fc m ,0) r (±2fc m ,0,0) 

3£>{(±2fc m ,±2fc m ,±2fc m ) 2d\ (±2k m ,0,±2k m ) ID i (0,±2ib m ,0) 

( 0 , ±2k m ,±2k m ) ( 0 , 0 ,± 2 * m ). 

( 31 ) 

By symmetry all the contributions within each class are equal. Therefore, 

£ (a,ias) (J k) = 6£S as) (ib) + 12 £S* s \k) + 8£$ as) (ifc) (32) 
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where £[ A p* s \k) is the contribution from any one of the ID modes, is 

the contribution from any one of the 2D modes, and ^ 3 ^ ias ^(A:) is the contribution 
from any one of the 3D modes respectively. If the modified wave-number, k, of 
a numerical method and the energy spectrum of the turbulence, E(k ), are known, 
(32) may be evaluated numerically using either E min (k) or F max (ib) to get the lower 
and upper estimates for the aliasing error respectively. 


2.2.3 Subgrid and total contributions 

The total nonlinear term N and the (exact) subgrid force S can be readily written 
down in terms of the Fourier-basis: 

Ni( k) = -iP,mn(k) j dk' dk"6(k' + k" - k)u m (k')«n(k"), (33) 

and 

5,(k) = -iP im „(k) (^j J - dk’dk"6(k' + k" -k)u m (k')a n (k"). (34) 

The power-spectra are defined as 

^F = v ! L m „T- ((S ' (k)S ' (k) ‘»= < 35 > 

^ { W< k ) iv '( k r))n (36) 

where { }n as usual denotes angular average over the sphere |k| = k. 

The evaluation of (36) is similar to the calculation of £ FD (Ar) in § 2.2.1. One only 
needs to replace ‘Ai mn ’ in (16) by i —Pi mn ' > and drop the last term involving the 
viscosity. The resulting expressions can be further simplified using the properties 
of the P % mn tensors (see paper 1): 

^(A:) = A ; 2 [F 1 (A:) + F 2 (A ; ) + F3(A:)]. (37) 


where Fi, F 2 and F 3 are as defined in (20). 

The computation of S(k) once again requires us to restrict the k space integration 
to a cubical domain which makes it difficult to handle the integrals analytically. This 
difficulty is dealt with in precisely the same manner as was done in the computation 
of the aliasing error. The cubical domain in k space is replaced by a spherical region 
of appropriate size. This is completely equivalent to replacing the energy spectrum 
E(k) by a pseudo-spectrum E min (k ) or F max (A:) defined as: 


^pmin 



E(k) if k > y/3 km 
0 otherwise 


(38) 


and 


c (jfc) = | if ^ > km 

1 0 otherwise. 


(39) 


With this modification, the calculation is exactly identical to that just presented 
for the nonlinear term. Thus, one obtains 


,S(fc) = A: 2 [F 1 (fc) + F 2 (A : )+F 3 (fc)]. (40) 

where in the evaluation of the functions F;, the pseudo-spectrum E mm (k ) or F max (fc) 
should be used in place of E(k) to obtain the lower and upper bounds respectively. 
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Wavenumber ( k ) 

FIGURE 1. The Von-Karman spectrum normalized so that the maximum energy 
density is at k = 1 and £(1) = 1. 


2.3 Application to LES 

The results established in the previous sections will now be applied to establish 
quantitative measures of errors in LES. In LES, the grid spacing A is typically much 
larger than the Kolmogorov length so that molecular viscosity plays a negligible role. 
Therefore V is set to zero throughout this section. For the energy spectrum we 
assume the “Von-Karman form” 


£ <‘> - iwr (41) 

where the constants a = 2.682 and b = 0.417 are chosen so that the maximum 
of E(k) occurs at k = 1 and the maximum value £*(1) = 1. This can always be 
ensured by a proper choice of length and time-scales. The Von-Karman spectrum 
has the property E{k ) ~ k 4 as k — ► 0 and E(k) ~ as k — > oo and is a fair 

representation of inertial range turbulence. A plot of this spectrum is shown in 

Fig. 1. 


2.3.1 Spectra 

The power spectra Af(k) and S(k) are evaluated numerically from (37) and (40), 
respectively, using the Von-Karman spectrum. We assume the LES filter to be 
equal to the grid spacing A. The results are shown in Fig. 2 for k m = 8 and 
32, where k m = 7r/A. It is seen that the power spectrum of the total nonlinear 
term is reasonably flat at high wavenumbers while the subgrid contribution rises 
monotonically to a maximum (which appears as a “cusp” when plotted on a linear 
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Wavenumber ( k ) 

FIGURE 2. The total nonlinear term ( ) compared to the lower ( A) and upper 

( ▼) bound of the true subgrid force for k m = 8 and 32. 


scale) at the cut-off wavenumber k m . The subgrid contribution is seen to be a 
relatively small part of the total contribution from the nonlinear term. 

Subgrid modeling is a very important part of large-eddy simulation. A parametriza- 
tion of the interaction of the unresolved eddies with the resolved ones is expressed 
sis a subgrid model. It is therefore desirable that the errors inherent in the nu- 
merical method be much smaller than the physically motivated subgrid model. We 
now examine to what extent such an expectation is realized for a second order 
central-difference method implemented with the nonlinear term in divergence form. 
A second order central-difference scheme is characterized by the modified wavenum- 
ber ki = sin(fc;A)/A (i = 1, 2 or 3). Eq. (23) is used to compute the power spectra 
of the finite-differencing error 8^ FD \k) for k m = 8 and 32. These results are com- 
pared to the power spectra of the respective subgrid terms in Fig. 3. Only two 
values of k m are shown for clarity. The figures have the same qualitative appear- 
ance for all values of k m . The power spectrum of the finite-differencing error rises 
to a maximum at k = k m in the same manner as the subgrid contribution. How- 
ever, for all values of k m the finite- differencing error is substantially larger than the 
subgrid contribution over the entire wavenumber range. 

Figure 3 indicates that the error in a second order scheme cannot be reduced to a 
level below the subgrid contribution by sufficiently refining the grid. As the grid is 
refined (A: m is increased), both the error as well as the subgrid force decrease for all 
wavenumbers. However, the error continues to dominate the subgrid force through- 
out the wavenumber range irrespective of the resolution. Let us now examine if this 
situation can be improved by using higher order central-difference schemes. Fig- 
ure 4 shows the finite-differencing error evaluated using (23) for a second, fourth, 
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FIGURE 3. Finite-differencing error for the second order central-difference scheme 
( ) compared to the lower ( A) and upper ( ▼) bounds of the true subgrid force 



Wavenumber (&) 

FIGURE 4. Finite-differencing errors ( ) compared to the lower (A) and upper 

(▼) bounds of the subgrid force for k m = 8. The numerical schemes considered 
are second (highest curve), fourth, sixth, and eighth (lowest curve) order central- 
differences. 
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FIGURE 5. The aliasing error for a second-order central-difference method ( ), 

a fourth-order central-difference method (- - -) and an undealiased pseudo- 

spectral method ( ), compared to the lower (▲) and upper (▼) bounds of the 

subgrid force. Each method is represented by a pair of curves corresponding to the 
lower and upper bounds of the error. 


sixth, and eighth order central-difference scheme together with the subgrid term, 
computed using (40) for a fixed resolution, k m = 8. It is seen that higher order 
schemes do lead to reduced levels of error. However, even with an eighth order 
scheme, the subgrid contribution is dominated by numerical errors in about half of 
the wavenumber range. 

Figure 5 shows the corresponding comparison for the aliasing error computed 
using (30). In general, increasing the order of a scheme has a relatively weak effect 
on the aliasing error and the effect is primarily in the high wavenumber region. This 
effect is in fact in the “reverse” direction compared to the finite-differencing error. 
That is, the lowest pair of curves which correspond to a second-order scheme have 
the smallest aliasing error and the highest pair corresponding to an undealiased 
pseudo-spectral method have the largest. The aliasing errors for sixth and eighth 
order schemes are intermediate between the fourth and the pseudo-spectral; they 
have been omitted from Fig. 5 for clarity. The effect is of course quite easy to 
understand. In the one dimensional problem, the aliasing part of the nonlinear 
term is multiplied by the modified wavenumber which approaches zero at the cut- 
off so that the aliasing error is also reduced to zero at k m . In the three dimensional 
problem a similar situation applies except that the power spectrum does not actually 
go to zero on account of the averaging over wavenumber shells. However, the aliasing 
error is reduced at high wavenumbers for central-difference schemes. 
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2.3.2 Scaling laws 

In this section, the dependence of some measure of “global error” on resolution, 
l m , is investigated. An appropriate measure of the kind is 


v, - 



(42) 


where V stands for TD\ ‘alias’, 4 nF or ‘sg’ corresponding to the global finite- 
differencing error, aliasing error, total nonlinear term, or subgrid term respec- 
tively. < 7 * is closely related but not exactly equal to the the rms value, which is 
given by the integral of the power spectrum over the entire wavenumber range. 
The correspondence is not exact because the modes at the corners of the cube 
k m \ x l—km,k m ] x [— k m ,k m ] outside of the inscribed sphere of radius km 
have not been included in the definition (42). Thus, <r* is a lower bound of the true 
rms value. The <7* can be evaluated els a function of k m by numerically integrating 
the power spectra £^*\k) presented earlier. 

Figure 6 shows the lower and upper bounds (measured by the corresponding <7*) 
for the subgrid force <r sg as a function of k m . The corresponding quantity for the 
total nonlinear term <r n i is also shown for comparison. The subgrid contributions 
are seen to obey a power law. A least square fit gives 

_ J 0.36 k^ 39 (Lower bound) , 4 o\ 

<7sg ( 0.62 A:^ 48 (Upper bound) 


The total nonlinear term also appears to follow a power law. A least square fit in 
this case gives 

<r n i = 1.04 C 97 . (44) 


The fitted curves (43) and (44) are shown in Fig. 6 as dashed and solid lines respec- 
tively. Thus, the relative subgrid contribution is (roughly) cr sg /<7 n i ~ A;” 0,5 , that 
is, the role of the subgrid model decreases at higher resolution. As an illustration, 
for an LES that resolves about a decade of scales beyond the energy peak, the rms 
value of the subgrid force, according to this formula, should be in the approximate 
range 11 — 19 % of the rms value of the total force. 

The following heuristic “derivation” (Tennekes & Lumley 1983) is sometimes 
given for the scaling of the subgrid term. The traceless subgrid stress is = 2 VfSij 
where v% is the eddy-viscosity and Su is the rate of strain. The rate of dissipation 
e = TijSij = Vt\S \ 2 is a constant according to the classical Kolmogorov argument. 
Therefore, fail ~ v t |S| y/evt. Now, it seems reasonable to postulate that v t is the 


product of the grid-spacing, A, and the rms velocity of the subgrid eddies, 
The latter can be estimated from the Kolmogorov spectrum 



^~\L 


* roc 

1/2 r yoo 

/ E(k)dk 
Jk m 

~ / k~ 5/3 dk 

Uk m 


1/2 


(km)- 1 ' 3 ~ A 1/3 . (45) 
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FIGURE 6. Global measure of the toted nonlinear term, <r„; ( • ) and subgrid force, 
a„ g ( lower bound: ▲ , upper bound: ▼) plotted as a function of the maximum 
resolved wavenumber, k m . The lines represent power law fits obtained by the least- 
square method. 

Thus, v t ~ AA 1 / 3 ~ A 4 / 3 so that |r tJ | ~ y/eui ~ A 2 / 3 . The subgrid force, which is 
the derivative of should then scale as |r ,j|/ A ~ A -1 / 3 ~ ( k m ) 3 ^ 3 . The scaling 
exponent (0.4-0.5) in (43) is reasonably close to what this rough argument predicts. 
It should be noted that, even though the subgrid stress decreases with increasing 
resolution, its derivative, the subgrid force, actually increases. 

Figure 7 shows the integrated finite-differencing error, <7 fdi plotted against k m . 
There appears to be an asymptotic approach to a power law behavior for large k m . 
A least square power law fit to the last three data points gives 

{ 1.03 (Order 2) 

0.82 (Order 4) 

0.70 (Order 6) (46) 

0.5 (Order 8) 

0 (Spectral) 

which are shown as solid lines in Fig. 7. The subgrid terms cr S g are also shown for 
comparison. It is significant that the exponent in the dependence of the integrated 
error on resolution in (46) turns out to be independent of the order of the scheme. A 
higher order scheme reduces the error only through a reduced prefactor multiplying 
the ~ fc^ 75 term. 

Figure 8 shows the integrated value of the aliasing error cr alias plotted against 
k m . The lines are power law fits to the data. Only the second order scheme and 
the pseudo-spectral scheme without dealiasing is shown. The curves for the fourth, 



Discretization errors in LES 


19 



k 


m 


FIGURE 7. Finite-differencing errors, <7 fd plotted as a function of the maximum 
resolved wavenumber k m (x) for central differencing schemes of order 2 (topmost), 
4, 6 and 8 (lowermost). The solid lines are least-square power law fits. Lower (A) 
and upper (T) bounds of the subgrid force cr sg are also shown for comparison. 



k 


m 


FIGURE 8. Upper and lower bounds of the aliasing error, cr a i; as , for a second-order 
(+) and undealiased pseudo-spectral method (x). Solid and dashed lines are least- 
square power law fits. Upper (▼) and lower (A) bounds for the subgrid term cr sg 
are also shown for comparison. 
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sixth, and eighth order schemes have intermediate positions and have been omitted 
for clarity. These fits are given by the following analytical expressions; 


^alias 


' 0.90 fc^ 46 
2.20 k°J 6 
0.46 k° m il 
k 1.29 k° m 65 


(Lower bound, Pseudo-spectral) 
(Upper bound, Pseudo-spectral) 
(Lower bound, Second-order) 
(Upper bound, Second-order) 


(47) 


The important distinction from Fig. 7 is that here the curves are “reversed”. Thus, 
the lowest curve corresponds to the second order scheme and the highest corresponds 
to an undealiased pseudo-spectral scheme. The subgrid term a sg is also shown for 
comparison. Of course, for a spectral scheme properly dealiased with the ‘3/2-rule’ 
both the aliasing as well as the finite-differencing errors are identically zero. 

2.3.S Discussions 

The results of the above analysis may be summarized as follows. In large eddy 
simulation, the net effect of the unresolved eddies on the resolved ones is repre- 
sented by a subgrid model. The resulting equations, which are the Navier-Stokes 
equations augmented by an additional term, the subgrid force, is then solved nu- 
merically. In such a procedure the presumption is that the associated numerical 
errors are small compared to the subgrid model being used. To keep the analy- 
sis as simple as possible, isotropic turbulence in a ‘box’ with periodic boundary 
conditions was considered together with a simple numerical method: an order n 
(n = 2 to 8) central-difference scheme with the nonlinear term in the divergence 
form. It was found that the power spectrum of the aliasing error is significantly 
larger than the subgrid term over most of the resolved wavenumber range. Higher 
order schemes have the effect of increasing the aliasing error. The finite- differencing 
error for a second-order scheme also remains significantly larger than the subgrid 
term over most of the resolved wavenumber range. The situation is improved by 
going to higher-order schemes. However, even for an eighth-order scheme, the error 
dominates the subgrid term for almost half of the resolved wavenumber range. An 
increase in grid resolution makes the errors increase faster than the subgrid force 
so that the situation cannot be improved by grid refinement as long as the cut-off 
is in the inertial range. 

We now consider a possible remedy for this difficulty. In LES the Navier-Stokes 
equations are first ‘filtered’ to remove all scales below some ‘filter-width’, A/. The 
resulting equations are then discretized on a grid of spacing A^. In order that the 
smallest resolved scales be representable on the grid, it is required that A g < A/. 
In practice one most often assumes A g — A/, to minimize computational cost and 
accepts the consequence that the “marginal” eddies may not be well resolved. As a 
matter of fact, this distinction between A g and A/ is often ignored and one speaks 
of ‘filter-width’ and ‘grid-spacing’ interchangeably. However, if one expects to ade- 
quately resolve all scales up to ‘A it is natural to require that ‘ A^’ be several times 
smaller than ‘A /’ (Rogallo &: Moin 1984). Thus, we are led to consider an LES with 
a filter- width A / performed on a numerical grid of spacing A g < A/. Clearly, in any 
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FIGURE 9. The finite-differencing error ( ) for a second order centred- difference 

scheme with A / = NA g for N = 1 (uppermost curve), 2, 4 and 8 (lowermost curve). 
The lower (A) and upper (▼) bounds of the subgrid force are shown for comparison. 
= 8 is held fixed. 
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FIGURE 10. The finite-differencing error ( ) for A / = 2A ? for a second 

(uppermost), fourth and eighth (lowermost) order central-difference scheme. The 
lower (A) and upper (▼) bounds of the subgrid force are shown for comparison. 
4 = 8 is held fixed. 
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such computation all Fourier- modes between &£ = 7r/A / and k 9 m = 7t/A ? must be 
held at very low amplitudes, for otherwise these “contaminated” modes would soon 
destroy the accuracy of computation of the modes (0, k ) through nonlinear inter- 
actions. This might be achieved naturally by the effective “dissipation range” of 
the eddy- viscosity. This may also be achieved by replacing the usual discretization 
of the Navier-Stokes equations by the following alternative (Lund 1995) 


dui 

dt 


6 _ . 6F[P] &F[r ti ] 6 S 

6 xfWiUj] Sx _ 6x . +u Sx k 6x k U " 


(48) 


where F[ ) represents a suitably designed filtering operation that reduces the am- 
plitudes of all modes in the range (fc£, k^) to zero or very small values. [Compact 
filters for finite-difference schemes that are close to a sharp low pass filter in Fourier 
space were first considered by Lele (1992). They have been used in the present 
context by Lund (Lund 1995 ).] The finite-differencing operator 6/Sxj is on the 
finer grid A^. The effect of this modification is easy to investigate in the present 
formalism. Thus, for a second order method, the c A’ in the expression for the mod- 
ified wavenumber need simply be replaced with A g . Figure 9 shows the result of 
such a computation for a second-order central-difference method with A g = Af/N 
where N = 1, 2, 4, and 8 for a fixed = 8. It is seen that with N = 8, the 
finite- differencing error is about one or two orders of magnitude below the subgrid 
term throughout the wavenumber range from k = 0 to k^ = 8. However, taking 
A^ = A//8 increases the number of grid points by a factor of 8 3 = 512 and the 
total computational cost (if the time-step, At is limited by the CFL condition so 
that At ~ A) by a factor of 8 4 = 4096. It may therefore be advisable to use instead 
a higher order scheme in conjunction with a grid that is finer than the filter-width. 
In Fig. 10 A g has been fixed at A//2 and the spectra of finite-differencing errors is 
plotted for a second, fourth, and eighth order scheme. It is seen that for an eighth 
order scheme the finite-differencing error is several orders of magnitude below the 
subgrid term. The increase in computational cost due to the refined grid is a factor 
of 2 4 = 16. Implementation of an eighth order scheme would also carry a penalty 
in terms of added cost. However, in view of the vastly increased accuracy, the addi- 
tional cost may be justified. In addition to reducing the finite-differencing error, the 
filtering scheme (48) completely removes the aliasing error. This is because modes 
k' and k" that ‘alias’ to a mode k must satisfy the relation k' + k" — k = a where 
a is a member of the “reciprocal lattice” A. Any component of the vector on the 
left of this equation can be at most k £ so that the left-hand side cannot exceed 
3fc£. Since at least one component on the right-hand side is 2 k^ or larger, the 
equation cannot be satisfied if 3 k^ < 2 that is, if A / > (3/2)A g there cannot 
be any aliasing errors. This is of course the well known “3/2 dealiasing rule” (see 
e.g. Canuto et al. 1988). 


3. Future plans 

The analysis presented in this report is kinematic in nature in the sense that 
the departure of the right-hand side of the Navier-Stokes operator from its ideal 
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value is investigated. The effect of this error on the dynamics of the solution and 
ultimately on the prediction of averaged quantities is unknown. However, in the 
light of the present findings that these errors are comparable in size to the subgrid 
term, a careful and systematic study is required before finite-difference methods 
can be considered reliable. Such a program of study should choose a specific flow 
for which reliable experimental data are available and for which issues such as 
sensitivity to initial and boundary conditions are reliably known to be unimportant. 
Numerical simulations should then be performed using both spectral and various 
finite-difference methods and the results compared to experiments and to each other. 
The effect of reducing errors using methods described in § 2.3.3 on relevant statistical 
averages should be studied. 

A study of this nature has recently been undertaken by Kravchenko and Moin 
(Kravchenko and Moin 1995). They used a channel flow spectral code that uses 
B-splines in the wall normal direction and trigonometric basis functions in the 
homogeneous directions. By replacing the wavenumbers by the modified wavenum- 
bers in the homogeneous directions they were able to mimic various finite difference 
schemes. Numerical experiments were run with various forms (divergence, rota- 
tional, skew-symmetric) of the nonlinear terms with staggered as well as nonstag- 
gered grids. Aliasing errors in general were found to have a very serious effect on 
the simulation causing the flow to laminarize in some cases, as might be expected 
in the light of the present analysis. The effect of aliasing errors on the simulation 
as well as their size was found to depend strongly on both the form of the nonlinear 
term as well as the order of the scheme. Aliasing errors had the most serious effect 
for (undealiased) pseudo-spectral methods, a result also consistent with the present 
study. The effect of aliasing errors on numerical simulations have also been studied 
by Blaisdell et al (1995), Zang (1991), Kim et ai (1987) and Moser et al (1982) 
among others using numerical simulations. 
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